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An  electrochemical  cell  is  a  multidisciplinary  system  which  involves  complex  chemical,  electrical,  and 
thermodynamical  processes.  The  primary  objective  of  this  paper  is  to  develop  a  linear  graph-theoretical 
modeling  for  the  dynamic  description  of  electrochemical  systems  through  the  representation  of  the 
system  topologies.  After  a  brief  introduction  to  the  topic  and  a  review  of  linear  graphs,  an  approach  to 
develop  linear  graphs  for  electrochemical  systems  using  a  circuitry  representation  is  discussed,  followed 
in  turn  by  the  use  of  the  branch  and  chord  transformation  techniques  to  generate  final  dynamic  equations 
governing  the  system.  As  an  example,  the  application  of  linear  graph  theory  to  modeling  a  nickel  metal 
hydride  (NiMH)  battery  will  be  presented.  Results  show  that  not  only  the  number  of  equations  are  reduced 
significantly,  but  also  the  linear  graph  model  simulates  faster  compared  to  the  original  lumped  parameter 
model.  The  approach  presented  in  this  paper  can  be  extended  to  modeling  complex  systems  such  as  an 
electric  or  hybrid  electric  vehicle  where  a  battery  pack  is  interconnected  with  other  components  in  many 
different  domains. 
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1.  Introduction 

Due  to  the  recent  interests  in  battery  electric  and  hybrid  elec¬ 
tric  vehicles,  a  significant  amount  of  research  has  been  focused  on 
secondary  batteries  or  electrochemical  energy  storage  devices.  For 
this  reason,  many  of  these  battery  works  have  been  developed  as  a 
part  of  simulation  models  of  these  vehicles.  These  works  are  some¬ 
times  based  on  empirical  relationships,  at  other  times  on  a  detailed 
description  of  the  physical  and  chemical  processes  that  take  place 
in  the  cell  [1-4],  and  even  on  the  development  of  equivalent  circuits 
[5,6].  Various  techniques  have  been  used  to  develop  these  models 
such  as  lookup  tables,  lumped  parameter  models  [4],  or  distributed 
models  using  porous  electrode  theory  [1-3]. 

In  this  paper,  we  propose  a  formalism  which,  we  believe,  is  more 
appropriate  for  the  phenomenological  description  of  electrochem¬ 
ical  systems  which  usually  consists  of  complex  phenomena  across 
multiple  domains;  namely  the  chemical  domain,  electrical  domain, 
thermal  domain,  and  other  domains  especially  when  the  battery  is 
placed  in  a  larger  system  such  as  a  hybrid  electric  vehicle  system. 
Modeling  engineers  usually  cope  with  the  generation  and  solution 
of  the  equations  governing  the  motion  of  such  systems. 

Linear  graph  theory  is  a  branch  of  mathematics  that  stud¬ 
ies  the  manipulation  of  topology  [7,8].  Although  this  theory  has 
been  extensively  incorporated  into  formulation  of  a  wide  range  of 
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physical  systems,  namely  electrical,  mechanical,  and  hydraulic  sys¬ 
tems,  the  extent  to  which  this  theory  has  been  applied  to  modeling 
electrochemical  and  thermal  processes  remains  from  nil  to  mini¬ 
mum.  It  is  the  goal  of  this  paper  to  examine  this  particular  problem 
in  some  detail.  It  will  be  shown  in  this  paper  that  the  electro¬ 
chemical  processes  and  thermodynamic  behaviors  of  batteries,  in 
general,  can  be  described  as  equivalent  electrical  components  inter¬ 
connected  to  each  other,  making  it  possible  to  use  graph  theory  to 
develop  the  dynamic  equations  for  the  whole  system. 

The  paper  begins  with  a  brief  overview  of  linear  graph  theory 
and  associated  mathematical  theorems,  followed  by  a  discussion  of 
the  applications  of  linear  graphs  to  modeling  electrochemical  cells. 
An  example  will  also  be  provided  to  demonstrate  the  use  of  linear 
graphs  to  model  a  NiMH  battery  including  the  thermal  effects  and 
side  reactions.  Finally  are  some  concluding  remarks. 

2.  Linear  graph  theory 

2.1.  Overview 

A  linear  graph  representation  of  a  physical  system  is  seen  as  a 
collection  of  oriented  line  segments  called  edges  which  intersect 
only  at  their  node  points.  Although  physical  systems  in  different 
energy  domains  use  different  interpretations  of  nodes  and  edges, 
the  linear  graph  topological  interpretations  of  these  systems  are  the 
same:  nodes  are  the  boundaries  of  a  component,  while  a  set  of  edges 
represent  the  component  itself.  For  example,  the  linear  graph  for 
the  electrical  network  given  in  Fig.  1  can  be  constructed  by  drawing 
a  node  for  each  point  at  which  two  physical  elements  connect,  and 
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Fig.  1.  Electrical  network  example. 


interconnections.  That  is,  the  system’s  topological  equations  are 
always  linear  and  may  be  formulated  in  a  systematic  fashion 
regardless  of  the  linearity/nonlinearity  of  a  system’s  component 
equations.  The  complete  topological  description  of  a  physical  sys¬ 
tem  can  be  written  in  a  simple  mathematical  form  using  an 
incidence  matrix  T.  This  is  a  v  x  e  matrix,  where  v  is  the  number 
of  nodes  in  the  graph,  and  e  is  the  number  of  edges.  The  incidence 
matrix  has  elements 

'  - 1 ,  if  edge  j  is  incident  upon  and  towards  node  i 
rr  _  J  +1>  if  edge  j  is  incident  upon  and  away  from  node  i 

^  —  |  0,  if  edge  j  is  not  incident  upon  node  i 

Specifically,  the  directed  edge  j  is  positively  incident  upon  node  i  if 
it  points  towards  the  node,  and  negatively  incident  if  it  is  directed 
away  from  the  node.  As  an  example,  the  incidence  matrix  for  the 
linear  graph  in  Fig.  2  takes  the  form: 


Fig.  2.  Linear  graph  isomorphic  to  electrical  network. 


by  replacing  these  elements  with  directed  edges  on  a  one-to-one 
basis.  This  linear  graph  is  shown  in  Fig.  2  and  is  said  to  be  topo¬ 
logically  equivalent  to  the  electrical  circuit  in  Fig.  1.  The  direction 
is  arbitrarily  assigned  to  each  edge  and  is  represented  by  an  arrow 
which  provides  a  reference  direction  for  the  two  abstract  variables 
associated  with  an  edge:  though  variable  and  across  variable. 

By  definition,  a  through  variable  is  a  variable  that  can  be  mea¬ 
sured  by  an  instrument  in  series  with  the  corresponding  element 
associated  with  the  edge,  while  an  across  variable  is  a  variable  that 
can  be  measured  by  an  instrument  placed  in  parallel  with  the  edge. 
For  an  electrical  network,  for  example,  through  and  across  variables 
can  be  the  currents  and  voltages,  respectively.  A  summary  of  possi¬ 
ble  through  and  across  variable  for  some  common  energy  domains 
are  summarized  in  Table  A.l .  These  two  quantities  are  carried  along 
the  entire  linear  graph  so  that  the  balance  of  energy  at  any  point  of 
the  graph  can  be  found.  This  makes  it  easier  to  deal  with  interfaces 
between  physical  systems  in  different  energy  domains  as  well  as 
determining  the  energy  within  the  system. 

The  through  and  across  variables  of  an  edge  are  not  indepen¬ 
dent  of  each  other,  but  are  related  by  a  mathematical  expression 
called  terminal  equation  or  constitutive  equation  which  represents 
the  physical  nature  of  the  linear  graph  component.  Clearly,  the 
empirical  nature  of  the  terminal  equation  associated  with  an  edge 
is  dependent  on  the  edge’s  domain.  For  example,  the  current  and 
voltage  associated  with  an  electrical  resistor  are  related  by  Ohm’s 
law  which  serves  as  the  terminal  equation  for  the  resistor  element. 

One  of  the  unique  features  of  linear  graph  theory  is  its 
ability  to  separate  the  equations  governing  the  physics  of  a  sys¬ 
tem’s  individual  components  from  the  equations  governing  their 


E\  i?2  E3  C4 
"  1  1  0  0  " 

0-111 
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in  which  the  vertices  associated  with  rows  and  the  edges  associated 
with  columns  have  been  explicitly  labeled. 

An  important  concept  in  linear  graph  theory  is  the  spanning  tree, 
which  is  a  subgraph  that  includes  all  the  nodes  of  the  original  graph 
without  any  loops.  The  remaining  edges  which  are  not  selected  in 
the  tree  are  grouped  in  cotrees.  The  edges  of  the  tree  are  called 
branches  while  the  edges  of  the  cotree  are  referred  to  as  chords.  In 
this  paper,  a  branch  is  represented  by  a  solid  line  while  a  chord  is 
depicted  by  a  dotted  line.  One  possible  tree  for  the  graph  in  Fig.  2 
has  been  drawn  with  solid  lines  and  consists  of  edges  E\  and  C4. 
The  remaining  edges  R2  and  R3  comprise  the  chords  of  the  cotree. 
For  graphs  consisting  of  multiple  parts  (as  in  systems  containing 
multiple  energy  domains),  each  part  is  independently  represented 
by  a  tree  and  the  collection  of  all  these  trees  make  up  the  graph’s 
forest,  while  all  of  the  chords  represent  its  coforest. 

Along  with  the  concept  of  a  system  tree,  two  new  topological 
matrices  can  now  be  introduced  -  the  fundamental  cutset  (f-cutset) 
matrix  and  the  fundamental  circuit  (f-circuit)  matrix.  A  cutset  is 
defined  as  a  set  of  edges  that,  when  removed,  divide  the  graph  into 
two  separate  parts.  An  f-cutset  consists  of  a  single  branch  and  a 
unique  set  of  chords.  On  the  other  hand,  a  circuit  is  a  set  of  edges 
that  form  a  closed  loop  with  the  f-circuit  being  a  circuit  containing 
one  chord  and  a  unique  set  of  branches.  In  an  electrical  system,  a 
cutset  is  essentially  a  linear  combination  of  the  node-based  Kirchoff 
current  law  (KCL),  which  is  an  expression  that  the  flows  passing 
through  a  node  are  conservative.  A  circuit  corresponds  to  the  Kir¬ 
choff  voltage  law  (KVL),  which  sums  up  the  forces  operating  along 
the  edges  enclosing  each  of  the  circuits.  Mathematical  representa¬ 
tions  of  the  f-cutset  (A f)  and  f-circuit  (B f)  matrices  can  be  written 
as 


Aft  = 

and 

BfCt  = 


1  Ac 

*b 

rc 

B  „  1 

CLb 

Ctc 

=  0 


=  0 


(1) 


(2) 


In  these  equations,  the  matrix  1  is  the  identity  matrix,  rb  and 
rc  represent  the  branch  and  chord  through  variables,  and  otb  and 
etc  corresponds  to  the  branch  and  chord  across  variables.  The  ele¬ 
ment  (i,  j)  of  the  matrix  Ac  takes  on  a  value  of  either  +1,  -1,  or  0 
which  indicates  whether  chord  j  is  a  part  of  and  oriented  in  the 
same  direction  as  the  defining  branch  i,  a  part  of  and  oriented  in 
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the  opposite  direction  of  the  branch  z ,  or  not  a  part  of  the  f-cutset, 
respectively.  In  a  similar  fashion,  the  value  of  the  element  (z,  j)  of 
the  matrix  Bb  is  either  +1,  -1,  or  0  depending  whether  branch  j  is  a 
part  of  and  oriented  in  the  same  direction  along  the  loop  as  chord 
z,  a  part  of  and  oriented  in  the  opposite  direction  along  the  loop  as 
chord  z,  or  not  a  part  of  the  f-circuit,  respectively. 

As  a  result  of  its  linearity,  the  cutset  equation  can  be  rearranged 
to  express  the  branch  through  variables  in  terms  of  the  chord 
through  variables.  This  arrangement  is  called  a  chord  transformation 
[7-1 1  ]  and  can  be  mathematically  expressed  as 

rb  =  -Actc  (3) 

In  a  similar  manner,  we  may  also  define  the  branch  transformation 
[7-1 1  ]  by  rearranging  the  circuit  equation  as 

clc  =  -Bbab  (4) 

Interestingly,  the  matrices  Ac  and  Bb  are  orthogonal  as  a  conse¬ 
quence  of  the  definition  of  f-cutset  and  f-circuit  matrices  [7,8,12] 

Ac  =  -B£  (5) 

2.2.  Formulation  of  system  equations  and  tree  selection 

The  f-cutset,  f-circuit,  and  incidence  matrices  can  be  used  to 
generate  the  governing  equations  for  the  physical  system  to  which 
the  linear  graph  is  topologically  equivalent.  In  the  branch-chord 
formulation  of  the  system  equations,  Eqs.  (3)  and  (4)  can  be  used 
to  eliminate  the  branch  through  and  chord  across  variables  from 
the  set  of  system  equations.  A  compact  set  of  final  equations  can 
be  obtained  by  substituting  (3)  and  (4)  into  the  terminal  equa¬ 
tions.  As  discussed  above,  the  cutset,  circuit,  and  terminal  equations 
provide  a  necessary  and  sufficient  set  of  equations  for  determin¬ 
ing  the  time  response  of  a  physical  system.  Thus,  the  selection 
of  trees  for  the  graphs  does  not  affect  the  underlying  mathemat¬ 
ical  model.  However,  the  selection  of  trees  can  greatly  reduce  the 
number  of  equations  that  have  to  be  solved  simultaneously,  espe¬ 
cially  if  some  care  is  taken  in  selecting  the  branches  of  these  trees. 
Leger  and  McPhee  [13]  made  an  observation  that  the  number  of 
dynamic  equations  remaining  will  depend  directly  on  the  number 
of  branch  coordinates  that  have  been  used.  Therefore,  the  num¬ 
ber  of  equations  can  be  reduced  further  by  selecting  into  the  trees 
those  elements  for  which  a  minimum  number  of  across  variables  are 
unknown.  The  result  will  be  a  smaller  number  of  branch  coordinates 
and  therefore,  a  smaller  number  of  final  equations.  This  observation 
can  also  be  extended  by  selecting  into  the  cotrees  the  edges  that  can 
minimize  the  number  of  chord  through  variables  so  that  the  num¬ 
ber  of  chord  coordinates  in  the  final  set  of  equations  is  reduced.  This 
approach  is  useful  when  we  want  through  variables  to  appear  in 
the  final  equations  and  will  be  demonstrated  in  an  example  given 
in  Section  4. 

Using  this  simple  criterion,  it  is  desirable  to  include  voltage 
sources  in  the  electrical  domain,  and  position  drivers  and  revolute 
joints  in  the  mechanical  domain,  into  the  trees  since  their  across 
variables  are  completely  known  functions  of  time.  Similarly,  cur¬ 
rent  sources  in  the  electrical  domain  and  force  actuators  in  the 
mechanical  domain  can  be  selected  into  the  cotrees  since  the  cor¬ 
responding  through  variables  appearing  in  the  dynamic  equations 
would  be  known  functions. 

3.  Linear  graph  models  for  electrochemical  cells 

Linear  graph  theory  has  been  extensively  applied  to  many 
physical  systems  in  different  energy  domains,  namely  mechani¬ 
cal  domain,  electrical  domain,  and  hydraulic  domain  [7-13].  Linear 
graphs  have  not  yet  been  used  to  describe  electrochemical  cells  and 
thermodynamic  processes. 


In  this  section,  a  graph-theoretic  representation  for  electro¬ 
chemical  systems  similar  to  the  circuit  diagram  in  electrical 
network  theory  will  be  introduced.  Aside  from  being  intuitively 
advantageous  in  presentation,  this  graphical  notation  will  reveal 
the  role  of  system  topology  in  dynamic  behavior.  We  will  present 
the  procedures  for  obtaining  the  dynamic  equations  governing  an 
electrochemical  cell  directly  from  the  graph,  and  consequently  one 
may  look  upon  the  network  graph  as  another  notation  for  the  dif¬ 
ferential  equations  themselves. 

3.1.  Batteries  and  electrochemical  processes 

Fig.  3  shows  a  schematic  of  a  typical  electrochemical  cell.  Every 
electrochemical  system  contains  two  electrodes  separated  by  an 
electrolyte  and  connected  via  an  electronic  conductor.  Ions  flow 
through  the  electrolyte  from  one  electrode  to  the  other,  and  the 
circuit  is  completed  by  electrons  flowing  through  the  external  con¬ 
ductor.  At  each  electrode,  an  electrochemical  reaction  is  occurring 
with  driving  forces  for  reaction  being  determined  by  the  thermo¬ 
dynamic  properties  of  the  electrodes  and  electrolyte.  In  general,  a 
chemical  reaction  on  an  electrode  can  be  written  as 

YskM?  ve~  (6) 

k 

where  Mk  is  the  symbol  for  the  chemical  formula  of  species  k,  sk 
is  the  stoichiometric  coefficient  for  species  k,  v  is  the  number  of 
electrons,  and  zk  is  the  original  charge  of  species  k.  For  example, 
consider  the  reaction 

Zn  +  20H-  ?=±  Zn0  +  2e“  +  H20  (7) 

In  the  above  chemical  equation,  sZn0  is  -1,  s0H-  is  2,  zZn0  is  0,  z0H- 
is  -1,  and  v  is  2. 

Following  historical  convention,  current  is  defined  as  the  flow  of 
positive  charge.  Thus,  electrons  move  in  the  direction  opposite  to 
that  of  the  convention  for  current  flow.  Current  density  is  the  flux 
of  charge,  z.e.,the  rate  of  flow  of  positive  charge  per  unit  area  per¬ 
pendicular  to  the  direction  of  flow.  The  behavior  of  electrochemical 
systems  is  determined  more  by  the  current  density  than  by  the  total 
current,  which  is  the  product  of  the  current  density  and  the  surface 
area  of  the  porous  electrode.  In  this  paper,  symbol  j  refers  to  current 
density. 
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Owing  to  the  historical  development  of  the  field  of  electrochem¬ 
istry,  we  use  over-potential  to  refer  to  the  magnitude  of  the  potential 
drop  caused  by  the  resistance  to  the  passage  of  current  and  open- 
circuit  potential  for  the  potential  between  two  battery  electrodes  at 
which  no  current  flows.  This  open-circuit  potential  <p(t)  is  derived 
from  the  Gibbs  free  energy  and  then  is  reduced  to  the  Nernst  equa¬ 
tion  as 


m  =  U  +  (T(t)-T0) 


du 

3  T 


(8) 


where  U  is  the  open-circuit  potential  at  standard  conditions  for  the 
electrode,  T0  is  the  reference  temperature,  (dU/d  T )  is  the  reversible 
heat  constant  for  the  reaction,  and  ck(t)  is  the  concentration  or 
molality  of  the  reactant  k.  The  molality  ck{t)  can  be  related  to  the 
electric  charge  qk{t)  using  Faraday’s  law: 


Cfc(t)  = 


g kit) 
2  F 


(9) 


while  the  relationship  between  charge  and  current  can  be  repre¬ 
sented  by 


ifc(t) = 


dqk{t) 

dt 


(10) 


The  other  parameters  F,  R,  and  T  are  the  Faraday  constant,  the  gas 
constant,  and  the  battery  temperature,  respectively. 

The  resistive  force  for  the  chemical  reactions  at  the  two  elec¬ 
trodes  is  termed  the  surface  over-potential  and  is  given  the  symbol 
p.  The  current  density  j(t),  which  is  directly  related  to  the  rate  of 
chemical  reaction,  can  be  expressed  as  the  function  of  the  surface 
over-potential  by  the  Butler-Volmer  equation  in  the  form 


j(t)  =  »'o 


exp 


x  exp 


=  2  i0  sinh 


(: ) 

Ea  (  1  1 

r  \m  To 

aF 


-  exp 


(~mm) 


)] 


(m)m) 


x  exp 


Fa  (  1  1  \ 

[R  \T(t)  T0J J 


(11) 


The  Arrhenius  equation  exp [(Fa/F)((l/F(t))- (1/T0))]  represents 
the  dependency  of  the  reaction  rate  on  the  battery  tempera¬ 
ture  T(t)  and  the  activation  energy  Ea.  A  positive  p (t)  produces 
a  positive  (anodic)  current.  The  derivation  and  application  of  the 
Butler-Volmer  equation,  and  its  limitations,  is  discussed  in  Chapter 
8  in  the  work  of  [3].  The  dimensionless  parameter  a ,  called  charge 
transfer  coefficient  is  an  additional  kinetic  parameter  that  relates 
how  an  applied  potential  favors  one  direction  of  reaction  over  the 
other.  It  usually  has  a  value  between  0.2  and  2.0.  The  parameter  z0 
is  called  the  exchange  current  density  and  is  analogous  to  the  rate 
constant  used  in  chemical  kinetics.  A  reaction  with  a  large  value  of 
i0  is  often  called  fast  or  reversible.  As  an  example,  the  relationship 
between  the  current  density j(t)  and  the  surface  over-potential  p(t) 
for  the  main  chemical  reaction  on  the  positive  electrode  in  a  NiMFI 
cell  is  graphed  in  Fig.  4. 

Besides  the  main  chemical  reaction  that  generates  most  of  the 
battery’s  current,  there  are  usually  several  side  reactions  happening 
in  the  cell  container.  The  empirical  nature  of  these  side  chemi¬ 
cal  reactions  depends  on  the  type  of  batteries.  For  example,  the 
side  reactions  are  the  hydrogen  evolution  and  absorbtion  reactions 
in  a  lead-acid  battery  and  the  oxygen  evolution  and  absorbtion 
reactions  in  a  NiMH  cell.  These  side  reactions  may  have  a  signif¬ 
icant  impact  on  the  battery  performance  and,  therefore,  modeling 
these  effects  is  desirable.  Mathematically,  modeling  side  reactions 
is  similar  to  that  of  the  main  reaction  since  they  are  all  chemical 
reactions. 


Fig.  4.  Dependence  of  current  density  on  surface  over-potential  on  positive  elec¬ 
trode  of  a  NiMH  cell  at  30  °C. 

(a)  C  (b)  R 

- 1  | -  — WA — 

Fig.  5.  Circuitry  representation  for  (a)  open-circuit  potential  equation  and  (b) 
Butler-Volmer’s  equation. 


3.2.  Electrical  circuit-based  representation  for  chemical  reactions 

To  understand  the  composition  of  the  model,  an  electrical  cir¬ 
cuit  with  equivalent  components  will  be  developed.  Through  the 
relations  between  electrical  components,  this  equivalent  circuit 
will  clearly  show  the  relationship  between  the  electrochemical 
equations  to  facilitate  the  application  of  linear  graph  theory  to  the 
system. 

The  first  equation  we  will  examine  is  the  open-circuit  potential 
equation  (i.e.,the  Nernst  equation),  Eq.  (8),  which  relates  the  elec¬ 
trode  open-circuit  potential  0(t)  to  the  molality  of  the  substances. 
Since  molality  and  electric  charge  are  directly  related  by  Faraday’s 
laws  in  (9),  Eq.  (8)  can  be  thought  of  as  a  nonlinear  electrical  capaci¬ 
tor  (Fig.  5a)  which  stores  free  energy,  integrates  the  current  density 
j(t)  in  order  to  obtain  the  electric  charge,  and  thereby  obtains  the 
open-circuit  potential.  Starting  off  from  the  initial  battery  state  of 
charge  which  is  defined  by  the  molality  of  materials  on  the  bat¬ 
tery  electrodes,  the  charge  gradually  increases/decreases  when  the 
cell  is  charged/discharged  and  so  does  the  stored  electrochemical 
energy.  The  chemical  kinetics  represented  by  the  Butler-Volmer 
equation  in  (1 1 )  can  be  modeled  by  a  nonlinear  electrical  resistor  in 
which  the  current  density j(t)  and  the  surface  over-potential  p(t)  are 
coupled  by  a  nonlinear  expression.  This  component  is  illustrated  in 
Fig.  5b. 

Since  the  current  density  in  Eq.  (11)  is  the  current  density  that 
moves  the  ions  (i.e.,the  molalities  in  Eq.  (8))  from  one  electrode  to 
another,  these  two  components  must  be  connected  in  series.  The 
total  electrical  potential  across  both  of  these  components  is  the  sum 
of  the  individual  potentials  across  each  of  the  two  components  and 
is  called  the  positive  or  negative  electrode  potential.  We  can  state 
that:  each  chemical  reaction  is  represented  by  one  pair  of  these  capac¬ 
itor  and  resistor.  If  we  only  consider  the  main  chemical  reactions, 
there  will  be  two  pairs  of  capacitor  and  resistor,  each  representing 
a  chemical  reaction  on  one  electrode. 
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Electrode 


Fig.  6.  Circuitry  representation  for  multiple  reactions. 


Now  let  us  consider  an  electrode  at  which  there  are  n  chemical 
reactions  (i.e.,both  main  and  side  reactions).  Since  these  reactions 
are  independent  of  each  other,  the  total  electrode  current  den¬ 
sity  JtotaiCO  is  obtained  by  adding  up  the  current  densities  of  the 
individual  reactions.  That  is 

Jtotal(0  —  Jl(0  +J*2(0  H - +jn(0  —  ACe11  •  (12) 

^surf 

where  Asurf  is  the  surface  area  of  the  porous  electrode  and  icell 
is  the  current  produced  by  the  battery.  The  voltage  across  each 
capacitor-resistor  pair  is  also  the  same 

0C1(O+^1(O  =  0C2(O  +  ^2(O  =  ■■■  =0C„(O  +  ^„(t)  (13) 

It  can  be  inferred  from  Eqs.  (12)  and  (13)  that  the  electrode  can 
be  represented  by  multiple  capacitor-resistor  pairs  hooked  up  in 
parallel  as  shown  in  Fig.  6. 

We  can  close  the  circuit  by  connecting  the  positive  and  negative 
terminals  to  a  current  source  or  an  external  load  as  shown  in  Fig.  7. 
The  external  load  could  be  a  complete  electric  vehicle.  In  this  figure, 
/?int  is  the  internal  resistance  of  the  battery.  For  some  batteries,  this 
resistance  is  very  small  and  its  effects  can  be  ignored. 

3.3.  Thermal  effects 

So  far,  we  have  presented  the  electrical  circuit  representation  for 
the  chemical  reactions  in  an  electrochemical  cell.  The  model  that 
we  consider  assumes  that  the  battery  temperature  is  constant  or, 
in  other  words,  the  battery  model  we  have  investigated  so  far  is  an 
isothermal  model.  This  assumption  is  generally  acceptable  for  small 
cells  where  the  applied  current  is  not  high.  However,  when  the  cur¬ 
rent  intensity  is  high,  as  in  the  case  of  traction  batteries  for  electric 
or  hybrid  electric  vehicles,  the  effects  of  battery  temperature  can 
become  significant. 

Application  of  energy  balance  [14,15]  to  the  whole  cell  yields 
Cpm  cell  ^  =  —  fiAcen(T(t)  —  Ta)  +  iCell(Ot,cell(f ) 

-Jjk(t)(<Pk(t)~T(t)^)  (14) 

k= 1  '  ' 

In  this  equation,  cp  is  the  heat  capacity  of  the  cell,  mcen  is  the  mass 
of  the  cell,  h  is  the  external  heat  transfer  coefficient,  Acell  is  the 
cell  container  external  surface  area,  Ta  is  the  ambient  temperature, 
and  n  is  the  number  of  chemical  reactions.  The  right-hand  side  of 


the  equation  consists  of  three  terms:  the  first  term  corresponds  to 
the  heat  exchange  with  the  outside  environment  through  the  cell 
container  walls  according  to  Newton’s  law  of  cooling,  the  second 
term  refers  to  the  irreversible  heat  arisen  from  ohmic  heating  for 
the  whole  cell,  and  the  last  term  is  the  reversible  entropic  heat 
released  or  absorbed  by  the  chemical  reactions.  Eq.  (14)  shows  that 
the  heat  generation  rate  is  equal  to  the  sum  of  heat  transferred  out 
of  the  system  and  the  heat  stored  in  the  system. 

According  to  classical  thermodynamics,  Eq.  (14)  can  be  written 
as 

dQ(t)  _  dQ.ext (t)  dQ[YY(t)  dQrev(t)  MS'! 

dt  ~  dt  +  dt  dt  ’  1  J 

in  which  dQ[t)  =  cpmcelldT(t)  according  to  the  definition  of  heat 
capacity  and  (dQext(f)/df)  =  -  hAceU(T(t)  -  Ta)  according  to  Newton’s 
law  of  cooling.  The  last  two  terms,  (dQ,irr(t) / dt)  =  ice \\{t)vceii{t)  and 
(dQrev(t)/dt)  =  -ELi mmt)  -  T(t)(dUk/dT)),  are  the  rates  of 
heat  dissipated/absorbed  due  to  the  internal  resistance  and  chem¬ 
ical  reactions.  In  electrochemistry,  Qrev(t)  is  also  called  the  Gibbs 
free-energy  change. 

In  order  to  develop  a  linear  graph  for  the  thermal  domain, 
we  need  to  transform  the  thermal  balance  equation  into  the 
temperature-entropy  form  so  that  the  through  S(t)  and  across  T(t) 
variables  appear  explicitly  in  the  equation.  Dividing  both  sides  of 
Eq.  (15)  by  T(t)  results  in 

S(t)  =  Sext(t)  +  SiTT(t)  +  Srev(t),  (16) 

where  Sext(0.  5irr(t),  and  Srev(f)  are  the  time  derivatives  of  the 
entropy  for  the  external  temperature  exchange  term,  the  irre¬ 
versible  term,  and  the  reversible  term,  respectively.  In  Eq.  (16), 
Srev(0  is  the  sum  of  the  individual  entropies  for  the  chemical  reac¬ 
tions  and  can  be  expressed  as  Srev(0  =  Ek=i^rev/<(f)* 


3.4.  Linear  graph  for  battery  model 

Following  the  circuitry  representation  for  the  chemical  domain 
as  discussed  in  Section  3.2  and  the  temperature-entropy  represen¬ 
tation  given  in  Section  3.3,  we  can  develop  the  linear  graphs  for 
both  domains.  Examples  for  such  graphs  are  shown  in  Fig.  8.  In  this 
figure,  C’s  and  fl’s  are  the  nonlinear  capacitors  and  resistors  whose 
equations  are  given  in  (8)  and  (11).  For  convenience,  the  current 
i(t)  and  voltage  v(t )  will  be  used  as  through  and  across  variables  for 
the  chemical  domain.  The  current  flowing  through  each  component 
can  be  related  to  the  current  density  by  z(t)=Asur£/(t).  The  voltage 
across  each  resistor  is  the  surface-over  potential  rj(t)  while  the  volt¬ 
age  across  the  each  capacitor  component  is  the  open-circuit  voltage 
0(t).  The  tree  branches  and  chords  have  been  arbitrarily  chosen  as 
shown  in  Fig.  8a,  which  bears  a  striking  resemblance  to  the  physical 
system  in  Fig.  7.  We  use  a  solid  line  to  represent  a  tree  branch  and  a 
dotted  line  for  a  chord.  If  the  model  of  the  external  circuit  is  known, 
we  can  also  construct  the  linear  graph  for  the  entire  system.  For  bat¬ 
tery  charge  and  discharge  operations,  the  external  circuit  is  simply 
a  current  source  which  delivers  electric  current  to  the  battery.  The 
two  graphs  in  Fig.  8  are  coupled  by  the  temperature  variable  T(t ) 
which  appears  in  Eqs.  (8)  and  (11). 

For  the  thermal  domain,  the  graph  is  simply  a  set  of  edges  con¬ 
nected  in  parallel.  The  through  and  across  variables  for  the  thermal 
domain  are  the  time-derivative  of  entropy  5(t)  and  the  battery  tem¬ 
perature  T(t),  respectively.  It  can  be  realized  that  the  product  of  5(  t) 
and  T(t)  is  power,  same  as  the  product  of  voltage  and  current.  This 
indicates  that  the  energy  flowing  through  the  system  components 
is  conserved. 
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Positive  electrode  Negative  electrode 


Fig.  8.  Linear  graph  representation  for  (a)  chemical  domain  and  (b)  thermal  domain. 


The  system  dynamic  equations  can  be  developed  following  the 
procedures  discussed  in  Section  2.2.  We  can  write  the  through  vari¬ 
ables  vector  for  the  chemical  domain  as 


icjt)  ...  iCan(t )  iRcl(t)  ...  iRan(t)  iRint(t)  ice„(t) 


T 


(17) 

The  molalities  ck's  in  Eq.  (8)  can  be  replaced  by  charge  vari¬ 
ables  qkf s  using  Faraday’s  law  in  (9).  These  charges  are  also  through 
variables 

r  i T 

Qcc,(t)  ...  qcan(t)  qR cl(t)  ...  qRan(t )  qRiJt)  qce„(r) 


(18) 

Similarly,  the  across  variable  vector  can  be  defined  as  follows 

r  T 

I 'Ccl(t)  •••  vCan(t)  vRcl(t )  ...  vRan(t )  t'Rinc(t)  Fceii(r) 

(19) 


We  can  also  define  the  through  and  across  variable  vectors  for 
the  thermal  domain 


s  = 


5ext(f)  5j  rr(t)  -Srev!  (t)  •  •  •  5revn  (0  S(t) 


(20) 


and 


Text(f)  T'irr  (0  ^rev^t)  •••  TreVn(t)  T(t) 


(21) 


The  branch  and  chord  transformations  can  be  applied  directly 
to  the  current  and  voltage  variables  as  shown  in  Eqs.  (4)  and  (3). 
However,  for  the  charge  variables,  initial  values,  which  appear  as 
we  integrate  the  current  variables,  must  be  included  in  the  chord 
transformation  equation.  This  gives 

[qh(t)  +  qft(0)]  =  -Aaj[qc(t)  +  qc(0)]  (22) 

The  formulation  procedures  for  an  electrochemical  system  can 
be  summarized  as  in  Fig.  9.  The  figure  depicts  the  steps  in  the  for¬ 
mulation  as  well  as  their  resulting  output  variables.  After  the  final 
step  we  obtain  2(n-  1)  +  1  ODEs  representing  the  charge-current 
relation  equations,  2 (n  - 1)  + 1  algebraic  equations  for  the  chemi¬ 
cal  domain,  and  1  ODE  for  the  thermal  domain.  It  can  be  observed 
that  2 (n  -  1)  + 1  is  also  the  number  of  chords  in  the  linear  graph. 
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Fig.  9.  Formulation  steps  for  electrochemical  cells. 


The  number  of  equations  can  be  reduced  further  if  some  of  the 
through  variables  are  known  functions  of  time.  As  an  example,  if 
the  applied  current  on  the  battery  terminals  is  known  then  icell 
can  be  considered  as  a  current  driver  and  two  equations  (z.e.,  one 
charge-current  relation  and  one  algebraic  equation)  can  be  elim¬ 
inated  from  the  final  equations  since  the  value  of  zcell  can  now  be 
substituted  directly  into  all  the  equations. 

4.  Application  to  NiMH  cell  model 


chemistry  together  with  side  reactions  and  thermal  effects  will  be 
presented,  followed  by  a  linear-graph-based  formulation  for  both 
chemical  and  thermal  domains. 

4.1.  NiMH  battery  chemistry 

The  chemical  reactions  on  the  two  electrodes  of  the  battery  can 
be  written  as  follows: 

Main  reaction  on  positive  electrode: 


discharge 

NiOOH  +  H20  +  e"  —  Ni(OH)2  +  OTT 

charge 

Side  reaction  on  positive  electrode: 
20H“  ->■  ^02  +  H20  +  2e“ 

Main  reaction  on  negative  electrode: 

discharge 

MH  +  OTT  ^  M  +  H20  +  e“ 

charge 

Side  reaction  on  negative  electrode: 


(23) 


(24) 


(25) 


—  02  +  H20  +  2e 


20H- 


(26) 


where  the  metal  M  in  the  negative  electrode  is  an  inter-metallic 
compound,  usually  a  rare  earth  compound.  During  charging,  oxy¬ 
gen  is  generated  at  the  nickel  electrode  and  the  gas  is  formed  when 
the  solubility  limit  in  the  electrolyte  is  reached.  The  oxygen  is  then 
transported  to  the  metal  hydride  electrode  where  it  is  reduced  by 
the  recombined  reaction  (26).  During  discharge,  the  oxygen  gen¬ 
eration  reaction  may  occurs  at  low  rates,  at  which  the  potential  of 
the  nickel  electrode  is  higher  than  the  equilibrium  potential  of  the 
oxygen  generation  [16,17]. 

The  electromotive  force  in  the  battery  as  defined  by  the  open- 
circuit  potentials  (z.e.,  Nernst’s  equations)  for  the  main  reactions 
(23)  and  (25)  on  the  positive  and  negative  electrodes  is 


01(0  =  l/i  +  (T(t)  -  T0)^l  +  In 


CH+,max  Qf+(0 
ceCH+(t) 


(27) 


03 (t)  =  u3  +{T{t)-T0)^2l  +  In  (<+  +  9.712x10- 


3T 


.2372 exp 


+0.2372  exp  (  --  lvimw  )  (28) 

y  QviH,max  J 

2.7302  x  10“4 

(CMH(t)/CMH,max)2  +  0.010768 

and  the  equilibrium  potential  for  the  oxygen  reactions  (24)  and  (26) 
is  given  by 


02(t)  =  u2  +  (T(t)  -  T0)^f  +  In 


(29) 


Due  to  its  dominance  in  almost  all  hybrid  vehicles  the  NiMH 
battery  model  has  been  chosen  to  demonstrate  the  technique  in 
this  paper.  The  NiMH  is  one  of  the  latest  battery  technologies 
and  has  many  advantages  over  the  other  more  commonly  used 
rechargeable  batteries  such  as  the  lead-acid  battery  or  the  nickel- 
cadmium  battery.  Some  of  these  advantages  include  higher  energy 
density,  more  environmental  friendly,  and  less  prone  to  memory 
(z.e.,  periodic  exercise  cycles  need  to  be  done  less  often).  The  NiMH 
battery  model  presented  is  a  modified  version  of  the  lumped  bat¬ 
tery  model  proposed  by  Wu  et  al.  [4].  In  this  section,  the  NiMH 


In  these  equations,  cH+(t)  is  the  concentration  of  Ni(OH)2,  cmh(0 
is  the  concentration  of  the  metal  hydride  (MH),  p02  is  the  par¬ 
tial  pressure  of  oxygen  gas,  and  T(t)  is  the  battery  temperature. 
Other  parameters  and  constants  are  listed  in  Table  B.l.  Eq.  (28)  was 
curve-fitted  from  the  experimental  data  of  a  nickel/KOH/LaNi5  bat¬ 
tery  using  a  nickel  oxide  positive  electrode  provided  by  Paxton  and 
Newman  [  1  ].  There  is  only  one  open-circuit  potential  equation  (29) 
for  the  two  side  reactions  (24)  and  (26)  since  side  reactions  are  cou¬ 
pled  together  by  the  oxygen  transport  from  the  positive  electrode 
to  the  negative  electrode. 
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Positive  electrode 


Negative  electrode 


The  rate  of  chemical  reactions  on  each  electrode  is  defined  by 
the  Butler-Volmer  equation  which  relates  the  current  density  jk(t) 
to  the  surface  over-potential  rjk{t)  by 


Jk(0  =  *o  ,k 


-  exp 


(- 


RT(t) 


(30) 


where  /<  =  1 ...  3  for  the  first  three  chemical  equations  (23)-(25). 
The  exchange  current  density  i0tk  for  each  chemical  equation  is 
given  by 


*0,1  —  *0,1, ref 

x  exp 

*0,2  =  *0,2, ref 
*0,3  =  *0,3, ref  ^ 

Fa  3 

x  exp  -M- 


cH+(f) 

°, 

x  CH+ ,  ref 

) 

Fa,  1  ( 

1 

R  \T(t) 

\  0.5  / 

— 

ce,ref  J 

V 

CmhW 

0, 

cMH,ref 

) 

(  1 

1 

~  7b 

0.5 


cH+,n 


cH+(f)\ 


0.5 


v^H+,max  ^H+,ref  / 


(31) 


Po2(t)\ 


P  o 


exp 


Ce 

t-e,  ref 


£«.2  /  1  1 

R  \T(t)  To 

0.5  /  ,  xX  0.5 


(32) 


CMH,max  ~  QytH(f)  \ 
0viH,max  —  0viH,ref  J 


(33) 


where  iotk,ref  is  the  exchange  current  density  at  a  reference  reactant 
concentration. 

For  the  oxygen  reduction  reaction  on  the  negative  electrode,  a 
limiting  current  equation  is  used  for  the  rate  of  reaction 


j4(0  =  - 


Po2(t), 


Po2,ref 

where 

*0,4  —  *0,4,  ref  ^p 


*0,4 


Fa,  4  I 

( 1 

Ml 

U(o 

To/J 

(34) 


(35) 


The  battery  current  icell(t)  can  be  calculated  from  the  charge 
balance  equations  on  the  electrodes  given  by 


*cell(0  —  >^post*pos^posC/l(0  +  J2(f))  (36) 

*cell(f)  —  — ^negQneg^negC/3(f)  +  J4(f))  (37) 


The  mass  balance  of  the  nickel  active  material  is  given  by 


Ji(t)  =  F 
J3(t)  =  F 


dcH+(f)  PNi(OH)2 

dt  PNi(OH)2^posQpos 

dcMn(t)  Tmh 

dt  PMH^neg^neg 


■4posOpos^posl2(f)  +  71negf*neg^negl4(f)  —  F  — 


dp0At)  y 


gas 


RT{t ) 


(38) 

(39) 

(40) 


The  battery  temperature  can  be  obtained  from  the  energy  bal¬ 
ance  of  the  whole  cell  described  by  Eq.  (14). 


4.2.  Linear  graph  representation  for  chemical  domain 

The  application  of  the  linear  graph  concept  to  the  chemical  reac¬ 
tions  of  the  NiMH  battery  is  a  straightforward  operation.  The  main 
and  side  chemical  reactions  for  the  NiMH  battery  model  shown  in 
Section  4.1  can  be  represented  by  an  equivalent  electrical  circuit 
as  shown  in  Fig.  10.  The  open-circuit  voltage  equations  (27),  (29) 
and  (28)  can  be  represented  by  nonlinear  electrical  capacitors  Q, 
C2,  and  C3.  The  relationship  between  the  electrical  potentials  </>/< 
and  concentrations  in  these  equations  is  similar  to  the  capacitive 
relationship  between  voltage  and  charge  in  an  electrical  capacitor. 

The  nonlinear  resistors  f?i,  R2,  and  R3  are  used  to  model  the 
resistive  relationship  between  the  current  density  and  overvoltage 
in  Eq.  (30).  The  current  density  for  the  oxygen  reduction  reaction 
in  Eq.  (34)  and  the  applied  current  at  the  battery  terminals  can 
be  represented  by  the  current  drivers  14(f)  =  I4  and  *ceii(f)  =  fapp.  For 
simplicity,  it  is  assumed  that  the  battery  is  cycled  with  a  constant 
current  Japp.  We  also  assume  that  the  battery  has  thin  electrodes 
and,  therefore,  we  can  ignore  the  influence  of  the  battery  internal 
resistance. 

The  linear  graph  for  the  chemical  domain  is  shown  in  Fig.  11. 
We  see  that  the  topological  graph  structure  comprises  eight  edges, 
to  which  we  have  assigned  arbitrary  sign  directions.  The  number  of 
equations  to  be  solved  simultaneously  can  be  reduced  by  selecting 
a  tree  and  using  a  branch-chord  formulation,  as  described  in  Section 
2.  Choosing  Ci,  C3,  C3,  R2,  and  R3  as  branches  and  f^,  /4,  and  /app  as 
chords  can  reduce  the  number  of  final  equations  as  it  will  be  shown 
that  the  equations  for  J4  and  /app  are  known  (i.e.,  limiting  current  for 
oxygen  reduction  reaction  and  constant  charge/discharge  current). 

The  column  matrix  of  through  variables  is: 
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Fig.  11.  Linear  graph  presentation  for  main  and  side  chemical  reactions. 


iCl(0  *C2(0  *C3(0  t )  43(0  4l(0  4(0  4ell(0 


(41) 

Vector  i  can  be  broken  down  into  the  branch  through  variable 
vector  ib  and  chord  through  variable  vector  ic  as 

T 

*Cl(0  *C2(0  *’C3(0  42(0  43(0 


4  = 


4l(0  4(0  4ell(0 


(42) 

(43) 


Since  the  currents,  which  are  the  derivatives  of  the  electrical 
charges,  are  through  variables,  the  charges  are  also  through  vari¬ 
ables  and  can  be  written  in  vector  forms  as 

n  T 

9ci(0  <to(  0  to(  0  QR2(t )  <ta(0 


Qb  = 


qc 


<4l(0  94(0  9cell(0 


l  T 


(44) 

(45) 


The  current  and  charge  variables  are  expressed  in  terms  of  the 
current  densities  and  concentrations  as  follows 

4l(0  —  ^pos0posOosJl(O 

42(0  —  ^pos  Opos  Iposj 2(0 

43(0  =  AnegQnegln  eg h(t)  (46) 

4(0  —  ^negQneg4egj4(0 


and 

9ci(0 

<ta(0 


E4PosLNi(OH)2 

PNi(OH)2 

FAnegLMH 


CH+(f) 


<?C2(0  +  Q4(t)  =  ^J^PoJP 


Pmh 

fVg 

RT(t)l'U2'' 


CMH(t) 


(47) 


Similarly,  the  across  variable  vector  can  be  defined  as  follows 


v= 


t'ci(t)  vC2 (t)  vC3(t)  VR2(t)  vR3(t)  vn(t)  v4(t)  vcen (t) 


vb  = 


Vci(f)  VC2(t)  Vc3(t)  VR2(t)  VR3(t) 


VRl(t)  V4(t)  Vcell(f) 


(49) 

(50) 


where  the  voltage  variables  are  related  to  the  battery  electrical 
potentials  by 


Vck(t )  =  0k(t) 

VRkU)  =  Vk(t)  k  =  1...3 


(51) 


For  convenience,  the  Butler-Volmer  equations  are  converted 
into  the  conductance  form  by  applying  the  inverse  hyperbolic  oper¬ 
ation  to  Eqs.  (30)  as 


RT(t ) 

^(f)  =  ^fTln 


m  . .  (Mt)\ 

2lo ,k  V  l  2i0  i 


+  1 


k  =  1 . . .  3 


(52) 


The  expressions  in  Eqs.  (27),  (29),  (28),  (34)  and  (52)  can  be 
written  as  functions  of  currents,  charges,  and  voltages  using  the 
relations  in  (46),  (47)  and  (51);  doing  so,  we  obtain  the  following 
set  of  six  terminal  equations 


4:i(0  = 


4:2(0  = 


Ui+im-To)^ 

|  OT(0  ln  /cH+,max 


U2  +  (T(t)-T0)^l 


EAposTNi(0H)2  -  PNi(OH)29ci(0^ 


cePNi(OH)29ci(0 

RT(t )  ^  f  y/((gC2(0  +  94(t))/(fVgas))^(0\ 


2F 


— 3 - ) 

VC3 (0  =  1/3  +  (T(t)  -  T0)^ +  ^  In  (ce2)+ 9.712x10 


+0.2372  exp 


,31/3  ,  *1(0, 

9r  f 

28.057pMH 

FAnegf'MH^MH,max 
2.7302  x  10“4 


to(0 


[((PMH/(f^neg^MHCMH,max))9C3(0]  +  0.010768 


4?2(0  = 


RT(t) 
0L2  F 


ln 


42(0 


(53) 


2i’o,2^pos  Qpos  Oos 


4*3  (0  = 


+  1 

RT{t ) 
0^3  F 


42(0 


2  4,2^pos  Opos  Oos 

43(0 


+  1 


In 


2io,3AiegQnegf] 


neg 


4*1  (0  = 


+  1 

RT(t) 
ot\ F 


43(0 


-V 

2io,3^negdIneg4eg  y 

4l(0 


+  1 


In 


24,l^posd[pos0os 


4i(0 


+  1 


2z'o,l  ^pos  Qpos  Oos 

In  a  similar  manner,  the  current  driver  in  (34)  can  also  be  written 


as 


WI— yyw*.. 

Po2,refCV  gas 


(54) 


(48) 
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We  now  define  the  f-cutset  and  f-circuit  matrices.  For  the  graph 
given  in  Fig.  1 1  and  the  given  tree  selection,  the  fundamental-cutset 
matrix  is  obtained: 


"1  0  0  0  0  -1  0  0  " 
0  1  0  0  0  1  0  -1 

0  0  1  0  0  0  1  1 

0  0  0  1  0  1  0  -1 

0  0  0  0  1  0  1  1 


(55) 


which  can  be  used  to  write  the  chord  transformation  equation  for 
this  system 


/ici(0\ 

iC2(0 

*C3(0 

42  (0 

V  43(f)/ 


"-1 

1 

0 

1 

0 


0  0  " 

0  -1 

1  1 

0  -1 

1  1 


/4i(0\ 

4(0 

\  4ell(0  / 


(56) 


Substituting  Eq.  (47)  and  the  first  equation  in  (59)  into  the  ter¬ 
minal  equations  in  (53)  we  obtain 


RT(t) 

aqF 


In 


4i(0 


2io,i^  posQposOos 


4i(f) 


2io,l^posQpos0os 


+  1 


+Ul+mt)_ro,^  +  ?TO 

/ FAposi-NifOHJjCH  .max  “  /°Ni(OH)2<7ci(t)\ 

\  Ce/0Ni(OH)2<?Cl(t)  )  + 


-U2+(T(t)-T0)-^ 


,  ,dU2  RT(t)  ( + 

-’iF  +  TT1" - Turn - 


RTjt) 

ol2F 


In 


42  (0 


2io,2^pos<3pos0os 


42  (0 


2io,2^posOposO°s 


+  1 


(60) 


=  0 


By  applying  the  chord  transformations  in  (56)  and  (57)  to  both 
current  and  charge  variables  in  (54)  and  (60),  we  obtain  the  follow¬ 
ing  set  of  two  equations: 


Aiegf*negWg[9cell(0  +  9c2(0)  ~  9fll(f)  +  9ftl(0)  —  9cell(0)  +  94(01^(0*0,4 


RT{t ) 


In 


FV  gasPo2,ref 


=  4(0 


4i(0 


2*0,l^pOsQpos0oS 


4i(0 


2*0,l^posQpos^] 


pos 


+  1 


+  1/1  +(T(0-7o)^-  + 


aqF 

RT(t)  / F71posFNi(OH)2CH+,max  “  PNi(OH)2  (9fll  (0  +  <4l  (0)  “  9ci(0))\ 

+  F1  Y  CePNi(OH)2(9i?i(0  +  9ki(0)  -  9ci(0))  y  + 

rr  t^2  RT(t)  f  y/mt){-qn{t)  +  9jn(0)  +  9cell(0  -  9cell(°)  +  9C2(0)  +  94(0) 

-U2  +  (i  (t)  -  lo)^-  +  «-> i—  ln  I  - « — 1  -  I  + 


(61) 


RTjt) 

0i2F 


ln 


3  T  '  2F 

-4l(0  +  *cell(0 

2*O,2^pos0posOos 


gas 


+  ' 


~ 4l(0  +  *cell(0 

2*0,2^pos<3pos^pos 


+  1 


=  0 


For  the  charge  variables,  initial  values  must  be  included  in  the 
chord  transformation  equations,  which  yields 

/ 9ci(0  -  9ci(0)\ 

9C2(0  -  9c2(0) 

9C3(0  -  9c3(0) 

QR2( 0  -  9^2(0) 

V  9i?3(0  -  9i?3(0)  / 

From  the  fundamental-cutset  matrix,  one  can  directly  obtain  the 
following  fundamental-circuit  matrix: 


0  0 

0  -1 

1  1 

0  -1 

1  1 


/  9ki(0  -  9ki (0)  \ 

94(0  -  94(0)  (57) 

\  9cell(0  —  9cell(0)  / 


It  should  be  noted  from  the  above  equation  that  icell(f)  is  the  cur¬ 
rent  applied  on  the  battery  terminals  and,  therefore,  is  completely 
known.  We  have  the  following  relationships 


4i(0 

4(0 


d9/n(0 

dt 

dq4(t) 

dt 


(62) 


*cell(0 


g9cell(0 

dt 


Therefore,  Eq.  (61 )  become  a  set  of  two  ordinary  differential  equa¬ 
tions  (ODEs)  which  only  consists  of  three  unknowns  q^i(t),  94(f), 
and  F(f)  and  can  easily  be  solved  using  a  numerical  integrator  with 
proper  initial  conditions,  if  the  temperature  is  known. 


Bf  = 


1  -1 
0  0 

0  1 


0 

-1 

-1 


-10  10  0 
0  -1010 

1  -1001 


(58) 


in  which  each  row  corresponds  to  an  equation  along  the  edges 
enclosing  each  circuit. 

Similarly,  making  use  of  Eq.  (4),  the  branch  transformation  for 
the  system  can  be  written  explicitly  as 


/  Vr\ 

v4 

\  ^cell 


) 


1 

0 

1 

0 

(yCl\ 

0  0-10-1 

VC2 

0  1-11-1 

VC3 

VR2 

- 

\Vr3  / 

(59) 


4.3.  Linear  graph  for  thermal  domain 

In  Section  4.2  we  have  shown  the  steps  to  develop  a  linear  graph 
and  system  equations  for  a  NiMFI  model  under  the  assumption 
of  constant  battery  temperature.  In  real  automotive  applications, 
studying  the  thermal  effects  of  batteries  is  of  particular  impor¬ 
tance  due  to  the  large  influence  of  battery  temperature  on  the 
battery  and  vehicle  performance.  Besides  reducing  battery  effi¬ 
ciency,  overheating  a  battery  may  even  cause  an  explosion  if  the 
battery  temperature  is  not  controlled.  Due  to  these  reasons,  it  is 
desirable  to  also  develop  a  linear  graph  for  the  thermal  effects  in  a 
car  battery. 

To  construct  a  linear  graph  for  the  thermal  domain,  we  can 
make  use  of  the  entropy-temperature  relationship  in  Eq.  (16)  for 
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Fig.  12.  Linear  graph  representation  for  thermal  domain. 


the  four  chemical  reactions.  Considering  the  time-derivative  of  the 
entropies  as  through  variables  and  the  temperature  as  an  across 
variable,  we  can  construct  the  linear  graph  for  the  thermal  domain 
as  shown  in  Fig.  12.  By  applying  the  branch  and  chord  transforma¬ 
tions  in  Eqs.  (59),  (56)  and  (57),  we  can  also  convert  Eq.  (16)  into  a 
function  that  is  only  dependent  on  T(t),  q^i(t),  and  q4(t). 

4.4.  Simulation  results 

The  model  developed  has  been  used  to  simulate  several  scenar¬ 
ios  in  order  to  observe  its  behaviors,  particularly  with  regard  to 
thermal  effects.  Table  B.l  contains  the  model  parameters  used  in 
the  simulations.  These  parameters  were  identified  using  homotopy 
optimization  from  a  3.4  A  h  NiMEI  battery  produced  by  North  Amer¬ 
ican  Battery  Company  (NABC)  based  on  the  reference  parameters 
provided  in  the  work  of  [4].  The  battery  data  was  measured  at  A&D 
Technology’s  laboratory  in  Ann  Arbor,  Michigan,  USA. 

The  battery  voltage  versus  time  for  four  different  discharge 
and  charge  rates  from  1  C  (3.4  A)  to  1/8C  (0.425  A)  are  shown  in 
Figs.  13  and  14  .  As  expected,  the  battery  voltage  drops/rises  more 
quickly  as  the  discharge/charge  current  is  increased.  Fig.  15  com¬ 
pares  the  battery  voltages  obtained  from  simulation  results  and 
battery  testing  data  at  a  constant  charge  and  discharge  rate  of  1  / 5  C. 
It  can  be  seen  that  these  results  are  in  good  agreement. 

The  battery  temperature  during  discharge  and  charge  versus 
time  is  shown  in  Figs.  16  and  17  .  It  is  assumed  that  the  battery 
has  been  cooled  down  to  room  temperature  (25  °C)  at  the  start  of 
each  simulation.  In  these  figures,  the  depletion  of  reactants  results 
in  high  over-potential  loss  that  causes  a  rapid  cell  temperature 
increase.  Due  to  the  temperature  exchange  with  the  environment, 
the  battery  temperature  is  flattened  out  at  the  end  of  the  cycle.  At  a 
high  charge  rate,  due  to  the  high  ohmic  and  over-potential  losses, 
high  charge  potential  is  needed  as  expected.  The  cell  heat  gener¬ 
ation  is  also  significant  at  high  discharge  currents  for  the  same 
reasons.  This  may  cause  the  cell  temperature  to  rise  more  than 
10°C.  At  high-rate  charge,  the  oxygen  generation  also  increases 
very  quickly,  contributing  significantly  to  the  increase  in  the  cell 
temperature.  Flowever,  during  discharge,  oxygen  gas  is  only  gen¬ 
erated  at  low  currents.  This  explains  the  difference  between  the  cell 
temperatures  in  the  two  figures  at  the  same  rates,  particularly  at  the 
high  rates.  Battery  temperature  control  is  therefore  very  important, 
particularly  in  a  battery  electric  vehicle  or  hybrid  electric  vehicle 
system  where  the  current  intensity  is  usually  very  high. 


Fig.  14.  Battery  charge  at  constant  rates. 


Time  Thrl 


| .  Experiment  Simulation  | 

Fig.  15.  Simulated  and  experimental  battery  voltages  at  1  /5  C  charge/discharge  rate. 
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. _ Time  [hr] _ 

| - 2  C - 1C .  1/2  C - 1/4  C| 


Fig.  16.  Battery  temperature  rising  from  initial  temperature  25  °C  at  constant  dis¬ 
charge  rates. 


Fig.  17.  Battery  temperature  rising  from  initial  temperature  25  °C  at  constant  charge 
rates. 

One  of  the  most  noticeable  results  is  that  the  linear  graph 
model  simulates  approximately  30%  faster  than  the  lumped  model 
(i.e.,all  equations  in  [4]  are  stacked  together)  at  all  charge  and  dis¬ 
charge  rates  as  summarized  in  Table  B.2.  The  simulation  times  were 
obtained  from  a  Dell™OptiPlex  2.9  GHz  desktop  computer  using 
Maple  dsolve  function  based  on  the  default  Runge-Kutta  Fehlberg 
numeric  integrator  with  the  same  settings  (abserr = relerr  =  10-7) 
for  both  linear  graph  and  lumped  models. 

5.  Conclusion 

In  this  paper,  we  have  presented  a  linear  graph  formulation 
for  systematically  generating  a  compact  set  of  dynamic  equations 
governing  electrochemical  systems.  By  carefully  managing  how 
equations  are  formed,  a  smaller  set  of  expressions  is  obtained. 
This  benefits  symbolic  implementation  by  reducing  the  size  of  the 
largest  expression  that  needs  to  be  handled  by  the  computer,  thus 
allowing  for  the  analysis  of  more  complicated  systems.  It  was  also 
shown  that  the  equations  obtained  using  linear  graph  theoreti¬ 
cal  approach  simulated  approximately  30%  faster  than  the  original 
lumped  model. 

Since  the  interconnections  between  a  system’s  components  are 
represented  by  a  linear  graph,  tree  selection  strategies  can  be  used 
to  determine  the  modeling  variables  for  the  system.  It  is  clear  that 


this  flexibility  can  provide  benefits  during  formulation  as  well  as 
simulation. 

This  approach  can  be  extended  to  modeling  a  more  complex 
system  such  as  a  battery  electric  vehicle  or  a  hybrid  electric  vehi¬ 
cle  within  which  a  battery  model  is  an  important  part.  This  is  a 
potential  for  future  research  since  modeling  individual  parts  of 
a  battery/hybrid  electric  vehicle  has  been  done  in  the  literature 
[9,12,18],  but  never  before  using  linear  graph  theory. 
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Appendix  A.  Through  and  across  variables  for  some 
physical  systems 

See  Table  A.l. 

Appendix  B.  NiMH  battery  parameters  used  in  simulations 

See  Tables  B.l  and  B.2  . 


Table  A.1 

Through  and  across  variables. 


Domain 

Through  (unit) 

Across  (unit) 

Electrical 

Current 

Voltage 

i  (A) 

v(V) 

Mechanical  translational 

Force 

Velocity 

F(N) 

v  (ms-1) 

Mechanical  rotational 

Torque 

Angular  velocity 

r  (Nm) 

cb  (rads-1) 

Hydraulic 

Volume  flow 

Pressure 

(f)v  (m3  s-1) 

p  (N  m-2) 

Thermal 

Entropy  flow 

Temperature 

dS/dt  (WK-1) 

T(K) 

Chemical 

Molar  flow 

Chemical  potential 

dN/dt  (mol  s-1) 

li  (J  mol-1 ) 

Table  B.l 

NiMH  battery  parameters. 

Parameter 

Unit 

Symbol 

Value 

Specific  electrode  area  of 
positive  electrode 

cm2  cm-3 

Qpos 

4000.0 

Specific  electrode  area  of 
negative  electrode 

cm2  cm-3 

Oneg 

3000.0 

Surface  area  of  positive 
electrode 

cm2 

Apos 

175.0 

Surface  area  of  negative 
electrode 

cm2 

Aneg 

100.0 

Thickness  of  positive  electrode 

cm 

Ipos 

3.3  xlO-2 

Thickness  of  negative  electrode 

cm 

Ineg 

2.8  xlO-2 

Loading  of  nickel  active 
material 

gem-2 

f-Ni(OH)2 

6.8  x  10-2 

Loading  of  metal  hydride 
material 

gem-2 

Tmh 

1.13xl0- 

Concentration  of  KOH 
electrolyte 

mol  cm-3 

Ce 

7.0  xlO-3 

Reference  concentration  of 

KOH  electrolyte 

mol  cm-3 

Q,ref 

1.0  x  10-3 

Maximum  concentration  of 
Ni(OH)2  in  nickel  active 
material 

mol  cm-3 

cH+,max 

3.7  x  10-2 

Reference  concentration  of 
Ni(OH)2  in  nickel  active 
material 

mol  cm-3 

cH+,ref 

0.5ch+  max 

Maximum  concentration  of 
hydrogen  in  metal  hydride 

mol  cm-3 

CiVIH.max 

1.0  X  10-1 

material 
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Table  B.l  ( Continued ) 


Parameter 

Unit 

Symbol 

Value 

Reference  concentration  of 
hydrogen  in  metal  hydride 
material 

mol  cm-3 

cMH,ref 

0.5CMH,max 

Reference  oxygen  pressure 

atm 

Po2,ref 

1.0 

Exchange  current  density  of 
reaction  at  reference 
reactant  concentration  for 
first  reaction 

A  cm-2 

*0,1, ref 

15.1  x  10“4 

Exchange  current  density  of 
reaction  at  reference 
reactant  concentration  for 
second  reaction 

A  cm-2 

*0,2, ref 

2.0  x  10“4 

Exchange  current  density  of 
reaction  at  reference 
reactant  concentration  for 
third  reaction 

A  cm-2 

*0,3, ref 

10.2  x  10“4 

Exchange  current  density  of 
reaction  at  reference 
reactant  concentration  for 
fourth  reaction 

A  cm-2 

*0,4, ref 

13.2  xlO-4 

Activation  energy  for  first 
reaction 

J  mol-1 

Ea,1 

10.0  x  103 

Activation  energy  for  second 
reaction 

J  mol-1 

Ea,2 

120.0  x  103 

Activation  energy  for  third 
reaction 

J  mol-1 

Ea,  3 

10.0  x  103 

Activation  energy  for  fourth 
reaction 

J  mol-1 

Ea,  4 

10.0  x  103 

Reversible  heat  for  first 
reaction 

VK-1 

dUi/dT 

-1.35  xlO”3 

Reversible  heat  for  second 
reaction 

VIC1 

dU2/dT 

-1.68  xlO”3 

Reversible  heat  for  third 
reaction 

VIC1 

dU3/dT 

-1.55  xlO”3 

Reversible  heat  for  fourth 
reaction 

VK-1 

dU4/dT 

-1.68  xl0“3 

Charge  transfer  coefficient 

OL\ 

0.5 

Charge  transfer  coefficient 

a2 

1.0 

Charge  transfer  coefficient 

Oi  3 

0.5 

Open-circuit  voltage 

V 

uu,  uu 

0.527,  0.458 

Open-circuit  voltage 

V 

u2 

0.4011 

Open-circuit  voltage 

V 

u3 

-0.8279 

Open-circuit  voltage 

V 

u4 

0.4011 

Gas  volume  in  NiMH  cell 

cm3 

Vgas 

1.0  x  10-1 

Density  of  nickel  active 
material 

gem-3 

PNi(OH)2 

3.4 

Density  of  metal  hydride 

gem-3 

Pmh 

7.47 

Reference  battery  temperature 

K 

To 

303.15 

Table  B.2 

Average  simulation  time  (in  s)  comparison  between  linear  graph  model  and  lumped 
model. 


Applied  current 

Linear  graph  model 

Lumped  model 

Discharge  1  C 

0.898 

1.115 

Discharge  1/2C 

0.922 

1.176 

Discharge  1  /4  C 

0.902 

1.121 

Discharge  1/8  C 

0.916 

1.198 

Charge  1  C 

0.904 

1.128 

Charge  1  \2  C 

0.919 

1.122 

Charge  1  /4  C 

0.918 

1.124 

Charge  1  /8  C 

0.920 

1.189 
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